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Abstract 

A density functional theory (DFT) of lattice fermion models is presented, which uses the single- 
particle density matrix jij as basic variable. A simple, explicit approximation to the interaction- 
energy functional W[y] of the Hubbard model is derived from exact dimer results, scaling properties 
of W[y] and known limits. Systematic tests on the one-dimensional chain show a remarkable 
agreement with the Bethe-Ansatz exact solution for all interaction regimes and band fillings. New 
results are obtained for the ground-state energy and charge-excitation gap in two dimensions. A 
successful description of strong electron correlations within DFT is achieved. 



PACS numbers: 71.10.-w, 71.15.Mb, 71.10.Fd 



I. INTRODUCTION 



First principles methods and many-body lattice models are the two main theoreti- 
cal approaches to the electronic properties of matter. From the first-principles perspec- 
tive, the major breakthrough in the last decades has been Hohenberg-Kohn-Sham's (HKS) 
density-functional theory (DFT) and the derived powerful methods of electronic-structure 
calculation.0 Despite their unparalleled success in an extremely wide variety of problems, 
current implementations of DFT have still serious difficulties in accounting for phenomena 
that involve strong electron correlations as observed, for example, in heavy-fermion mate- 
rials, Mott insulators or high-T^ superconductors.! Being in principle an exact theory, the 
limitations of DFT have to be ascribed to the approximations used for the interaction-energy 
functional [p(r)] and not to the underlying formalism. The development of new functional 
improving the description of strong correlation effects is therefore a major current theoretical 
challenge. 

On the other side, the physics of strongly-correlated Fermi systems is intensively studied 
in the framework of parametrized lattice models (e.g., Hubbard, Anderson, etc.) by using 
specific leading-edge many-body techniques.! Taking into account the universality of DFT, 
and its demonstrated efficiency in complex ab initio calculations, it is quite remarkable 
that only few investigations have been concerned so far with applying the concepts of DFT 
to the lattice models describing strongly correlated fermions.lii'iil In fact, already from a 
formal standpoint, one may expect that DFT with an appropriate Ansatz for W should be 
a particularly valuable many-body approach to lattice models, thus becoming a subject of 
theoretical interest on its own. Moreover, DFT studies on simpler universal models also 
provide useful new insights relevant to first principles calculationsjil particularly since in 
some cases the exact solution of the many-body problem is available.! 

The purpose of this paper is to extend the scope of DFT to the description of strong 
electron correlations in lattice Hamiltonians and to demonstrate quantitatively for the first 
time the performance of lattice density-functional theory (LDFT) in one-dimensional (ID) 
and two-dimensional (2D) systems. Sec. II presents concisely the basic formalism of LDFT. 
In this framework the ground-state properties are obtained from the solution of exact self- 
consistent equations that involve derivatives of the interaction-energy functional W['~f] with 
respect to the single-particle density matrix 7. In Sec. |ITT| the dependence of W on the 
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nearest-neighbor (NN) density-matrix element 712 is analyzed and a simple explicit ap- 
proximation to 1^(712) is derived for the Hubbard model. Sec. |rV| discusses representative 
applications of this Ansatz. First, the accuracy of the method is demonstrated by compari- 
son with available exact results on the ID Hubbard model. New results are then discussed, 
particularly concerning the ground-state energy and charge-excitation gap in 2D lattices. 
Finally, Sec. |V] summarizes our conclusions. 



II. LATTICE DENSITY-FUNCTIONAL THEORY 

In order to be explicit we focus on the Hubbard model which is expected to capture the 
main physics of lattice fermions in a narrow energy band. The Hamiltonian 

H = tijclcja + f^il'^^h (1) 

includes nearest neighbor (NN) hoppings tij, and an on-site interactions U {ni„ = c\^Cia). 
The importance of electron correlations is controlled by one parameter, namely, the ratio 
U / 1. The hopping integrals Uj are defined by the lattice structure (typically, tij = — t < for 
NN ij) and thus play the role given in conventional DFT to the external potential Vext(^)- 
Consequently, in LDFT the single-particle density matrix 7jj replaces the density p(r) as 
basic variable, since the hopping integrals are nonlocal in the sites.S The situation is 
similar to the density-matrix functional theory proposed by Gilbert for the study of nonlocal 
pseudopotentials.lBi 

The ground-state energy Egg and density-matrix '-yf^ are determined by minimizing the 
energy functional 

E[^]=Ek[i] + W[^] (2) 

with respect to 7jj.@ The first term 

EK[l]=Y.'^ijlij (3) 

is the kinetic energy associated with the electronic motion in the lattice. The second term 
is Levy's interaction-energy functional^ given by 



WYf\ = niin 



(4) 
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where the minimization runs over all A^-particle states |^I/[7]) satisfying 

(^[7]lE4c,.|v&[7])=7., (5) 

for all i and j.&ii'i W[y] represents the minimum interaction energy compatible with a given 
7ij. It is a universal functional in the sense that it is independent of tij. However, note that 
it still depends on the type of model interaction, on the number of electrons A^e or band 
filling n = N^/Na, and on the number of sites iV^.E 

E['y] is minimized by expressing ■jij = 7ijt+7ijj, in terms of its eigenvalues rjka (occupation 
numbers) and eigenvectors Uika (natural orbitals) as 

k 

Introducing Lagrange multipliers fi and Xka {^ka = ^ka/Vka) to impose the usual constraints 
J2kcrVka = and J2i l^ifco-P = 1, one obtains the eigenvalue equations 

^ f dW\ _ 

/ , I tij + — I Ujka — ^ktjUika ( ' j 

with < /i {^ka > /i) if Vka = 1 iVka = 0), aud 

efc„ = /i if < ?7fc<^ < 1 . (8) 

In Eq. (^ self-consistency is implied by the dependence of dW/d'-yija on rjka and Uika- The 
present formulation is analogous to well-known results of density-matrix functional theory 
in the continuum.l However, notice the fundamental differences with respect to the KS-like 
approach proposed in Ref.@, which assumes non-interacting f-representability, and where 
only integer occupations are allowed. The importance of fractional orbital occupations to the 
description of electron correlations within density-matrix functional theory has already been 
stressed by Gilbert .0 In particular for the Hubbard model, one observes that < rjka < 1 for 
all k, except in very special situations such as f//t = or the fully-polarized ferromagnetic 
state. This can be understood from perturbation-theory arguments — none of the r/fco- is a 
good quantum number for U/t ^ — and has been explicitly demonstrated in exact solutions 
for finite clusters or the ID chain.0 Therefore, the case (|^) is the only relevant one in general 
and all Eka in Eq- (0) must be degenerate. Consequently, 

tij + = (9) 



for all i and j. Note that approximations of W in terms of diagonal 'ju alone can never yield 



At this point it is important to observe that the general functional, valid for all lattice 
structures and for all types of hybridizations, can be simplified at the expense of universality 
if the hopping integrals are short ranged. For example, if only NN hoppings are considered, 
Ek is independent of 'jij for pairs of sites ij that are not NN's. In this case, the constraints 
(^[7] I Z^o- cJo-Cjo-l^l'H]) = lij in Eqs. (H) and (||) need to be imposed only for i = j and for 
NN ij. This allows to reduce drastically the number of variables and simplifies considerably 
the search for practical approximations to W. Moreover, in periodic lattices the ground- 
state 'yfj is a translational invariant. In order to determine Egg and 'jfj, one may then set 
'Ju = n = Ne/Na for all sites i, and jij = 712 for all NN pairs ij. Thus, the interaction 
energy can be regarded as a simple function 1^(712) of the density-matrix element between 
NN's. It should be however noted that this also implies that W loses its universal character, 
since the NN map and the resulting dependence of W on 712 are in principle different for 
different lattice structures.! 

III. INTERACTION-ENERGY FUNCTIONAL FOR THE HUBBARD MODEL 

Given a self-consistent scheme that implements the variational principle, the challenge 
is to find a good, explicit approximations to the interaction-energy functional. ^^[7] may 
be determined exactly for small clusters by using numerical methods that perform the con- 
strained minimization explicitly.! For a Hubbard dimer with N^. = Na = 2 a. straightforward 
analytical calculations yields 



which represents the minimum average number of double occupations for a given degree of 
electron delocalization, i.e., for a given 712 (f/ > 0). Despite its simplicity, Eq. (|T^) already 
includes the fundamental interplay between electron delocalization and charge fluctuations, 
and provides useful insights on several general properties of ^^(712) that are valid for arbi- 
trary lattices: 

(i) The domain of definition of 1^(712) is limited by the pure-state representability of 712. In 
fact, 7i2 < 7i2 = 1; where 75^2 corresponds to the extreme of the kinetic energy (maximum 



such a behavior.! 




(10) 
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degree of delocalization) and thus to the U = ground-state for a given lattice and a given 
n. 

(ii) For 7i2 = 7i2) the underlying electronic state ^l/[7i2] is a single Slater determinant 
and therefore ^^(712) = -E'HF = n'^U/A, where -Ehf is the Hartree-Fock energy. Moreover, 
dW/d'-yu = 00 for 712 = 712, since 7^2 < 7i2 already for arbitrary small U/t, as expected 
from perturbation theory. 

(iii) Starting from 712 = 712, ^^(712) decreases monotonically with decreasing 712 reaching 
its lowest possible value, W = 0, for 712 = 7^ (7^ = for n = 1). The fact that W 
decreases with decreasing I712I shows that the correlation-induced reduction of the Coulomb 
energy is obtained at the expense of electron delocalization. 

(iv) 7^ represents the largest NN bond order that can be obtained under the constraint 
of vanishing Coulomb energy. A lower bound for 7^ is given by the bond order 7f2'^ in 
the fully-polarized ferromagnetic state which is formed by occupying the A^e lowest single- 
particle states of the same spin (n < 1). Note that the ground-state 712 always satisfies 
I12 ^ 7i2 < 7i2 even though, for n 7^ 1, it is possible to construct Ag-electron states having 



In order to derive a simple approximation to W^(7i2) that preserves the previous general 



that 1^(712) depends weakly on N^, Na and lattice structure if it is measured in units of E^p 
and if 712 is scaled within the relevant domain of representability [7]^, 712]- Physically, this 
means that the relative change in W associated to a given change in the degree of electron 
localization = (712— 7S)/(7i2~7i2) can be regarded as nearly independent of the system 
under study. A good general approximation to 1^(712) can then be obtained by applying 
such a scaling to the functional dependence extracted from a simple reference system which 
already contains the fundamental relationship between localization and correlation. We 
therefore derive an approximate ^^(712) taking its functional dependence from the exact 
result for the Hubbard dimer given by Eq. (0). In this way one obtains 



where -Erf, 7i2 and 7J5 are system specific [see (i)-(iv) above]. In practice, 7^^ may be 
approximated by the ferromagnetic fully-polarized 7f2'^ which is calculated, as 7^2, by inte- 
gration of the single-particle spectrum. 



712I < |7i~|. 



properties we 



take advantage of its scaling properties. Exact numerical studied^ have shown 




(11) 
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FIG. 1: Interaction energy ^^(712) of the one-dimensional (ID) Hubbard model as a function of 
the degree of electron delocalization (712 — 712)/ (712 ~ 7i2)- The symbols refer to exact results for 
different band fillings n and the solid curve to Eq. ([Tl|). 

Fig. |T] compares Eq. (|TT|) with the exact Wex{li2) of the ID Hubbard chain which is 
derived from the Bethe-Ansatz solution.0 One observes that the proposed approximation 
follows Wex{'^i2) quite closely all along the crossover from weak correlations (large W/U and 
712) to strong correlations (small W/U and 712)- This is remarkable, taking into account 
the simplicity of Eq. (0) and the strong band-filling dependence of -Ehf, 7i2) 7^. 
The quantitative discrepancies between Eq. (|ll]) and W(.x{li2) remain small in the complete 
domain of representability of 7 and for all band fillings: \W — Wexl/Eup < 0.063 for all 712 
and n. Consequently, a good general performance of the method can be expected already 
at this stage. In the following section several applications of LDFT are discussed by using 
Eq. (|lT]) as approximation to the interaction-energy functional. 

IV. RESULTS AND DISCUSSION 

In Fig. ^ the ground-state energy Egg of the ID Hubbard model is given as a function 
of band filling n for different values of Coulomb repulsion U/t. Comparison between LDFT 
and the Bethe-Ansatz exact solution shows a very good agreement. It is interesting to 
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FIG. 2: Ground-state energy Egg of the ID Hubbard model as a function of band filling n for 
different Coulomb repulsions U/t. The solid curves refer to the present lattice density-functional 
theory (LDFT) and the symbols to the Bethe-Ansatz exact solution.0 

observe that the accuracy of the calculated Egg is not the result of a strong compensation of 
errors since a similar accuracy is achieved for the kinetic and Coulomb energies separately. 
Indeed, as shown in Fig. ^, both local moments Sf = 3{{hii — hiiY) and kinetic-energy 
renormalizations are also very well reproduced as a function oiU/t. Moreover, notice that no 
artificial symmetry breaking is required in order to describe correctly the correlation-induced 
localization, as it is often the case in other approaches (e.g., antiferromagnetic spin-density 
wave for n = 1). For n < 0.8, the LDFT results are almost indistinguishable from the exact 
ones. Even the largest quantitative discrepancies, found for n = 1 and intermediate U/t, 
are acceptably small (e.g., \Egs — Eg^\/t = 0.044 for U/t = 4). For U t and ri = 1 we 
obtain Egg ~ —at^/U with a ^ 3.24 while the exact result is a = 41n2 ^ 2.77. The error 
in the coefficient a can be corrected by including in Eq. (pH]) a 4th-order term in gi2 which 
provides in addition with a systematic improvement for all values of the interaction strength 
{\El^, - Egg\/\E'g^^\ < 0.02 for all U/t)E 
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FIG. 3: (a) Kinetic energy Ek and (b) local magnetic moments Sl = ^{{ni-\ — niif') of the ID 
Hubbard model as a function of Coulomb repulsion U jt for different band fillings n as indicated in 
the inset of subfigure (b). The solid curves correspond to the present LDFT and the symbols to 
the Bethe-Ansatz exact solution.0 

Fig. ^ shows Egs of the 2D square lattice as a function of \J jt for representative band- 
fillings n. The LDFT results cover the complete range of model parameters involving es- 
sentially analytic calculations. As shown in Fig. good agreement is obtained with far 
more demanding ground-state quantum Monte Carlo (QMC) studiesEl for U jt = 4. The 
reliability of LDFT in 2D systems is confirmed by comparison with exact Lanczos diago- 
nalizations on small clusters of the square and triangular lattices. In the inset of Fig. |^ 
we consider for example a A^^^ = 3 x 4 cluster of the square lattice with periodic bound- 
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U/(U+4t) 

FIG. 4: Ground-state energy Egg of the Hubbard model on the 2D square lattice as a function 
of Coulomb-repulsion strength U/t and for different band fillings n. The solid curves refer to the 
present LDFT and the crosses to quantum Monte Carlo (QMC) calculations .0 

ary conditions and Ng = Na. Like in the ID case, the overall performance is very good, 
with the largest quantitative discrepancies being observed for intermediate values of U/t. 
For instance, for U/t = 1 one obtains \Egs — Eg^\/\Eg^\ = 4.4 x 10~^, and for U/t = 4 
\Egs — Eg^\/\Egg \ = 9.8 X 10^^. Results with similar precision are found for a the triangular 
2D lattice. In this case, using also a A^^^ = A^e = 3 x 4 cluster with periodic boundary condi- 
tions, we find \Egs- E^^^l/lE^^^] = 1.7 x 10"^ for U/t = 1, and - = 6.6 x 10^^ 
for U/t = 4. For both lattice structures \Egs — E^^l decreases quite rapidly away from 
half-band filling as in the ID chain (see Fig. H). LDFT, combined with Eq. ( pT]) for W{'yi2), 
provides a correct description of electron correlations in different dimensions and lattice 
structures. 

The charge-excitation or band gap 

AE, = Eg,{N, + 1) + Egs{N, - 1) - 2EgsiN,) (12) 

is a property of considerable interest in strongly correlated systems which measures the 
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FIG. 5: Ground-state energy Egs of the Hubbard model on the 2D square lattice as a function of 
band filling n for U/t = 4. The solid curve refers to the present lattice density- functional theory 
(LDFT) and the crosses to quantum Monte Carlo (QMC) calculationsEl. In the inset figure LDFT 
is compared to exact Lanczos diagonalizations for a Na = 3x4 cluster of the 2D square lattice 
with periodic boundary conditions. Results are here given as a function of U /t at half-band filling 
(n = N,/Na = 1). 



insulating or metallic character of the electronic spectrum as a function of U/t and n. It can 
be directly related to the discontinuity in the derivative of the kinetic and correlation energies 
per site with respect to electron density ?t,.0 Therefore, the determination of AEc constitutes 
a much more serious challenge than the calculation of Egg, particularly in the framework 
of a density-functional formalism. At half-band filling, AEc = in the uncorrelated limit 
{U/t = 0) and it increases with increasing U/t. For fZ/t — > oo, AEc ^ U + E^ where E^ is 
the energy of the bottom of the single-particle band {Ei, = — 4t for a ID chain and Eb = —8t 
for the 2D square lattice). Fig. ^ presents LDFT results for AEc in ID and 2D Hubbard 
models (n = 1). Comparison with the Bethe-Ansatz result Ji3 and with available QMC 
calculationsEl shows a good overall agreement. However, a more detailed analysis reveals 
that in the ID case the gap is significantly overestimated for U/t <^ 1. Here we obtain 
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FIG. 6: Charge excitation gap AEc of the 2D Hubbard model (square lattice, n = 1). In the inset 
results are given for the ID chain. The solid curves refer to the present LDFT and the crosses to 



AEc oc {U/tY, while the exact solution shows that for the infinite chain AE^ increases 
much more slowly, namely, exponentially in —t/U. This discrepancy reflects the difficulty 
to describe long range-effects using an interaction-energy which functional dependence is 
derived from the dimer. Thus, it is possible that a similar overestimation of the gap at small 
U /t may also affect our results on 2D lattices. For larger U /t the accuracy improves rapidly 
as electron localization starts to set in, and the relative error in AE^. vanishes. Therefore, 
the development of a Mott insulator with increasing U/t is described correctly. 

Finally, we would like to comment briefly on a few other applications: (i) Dimerized chains 
with hoppings t ± 6t have been investigated by allowing for alternations of 712 in Eq. ([TT]). 
One observes that the precision of the results improves systematically with increasing dimer- 
ization. For example, for Na = Ne = 12 and U/t = 4, we find \Egs- E^^\/\E^^\ = 3.3 x 10"^ 
1.4 X 10^^, and 2.6 x lO^'^ for 6t/t = 0, 1/2, and 3/4, respectively. The non-dimerized case, 
shown in detail in Fig. H, is in fact the most difficult one, since for a collection of dimers 



QMC calculations (2D, U/t = 4) or to exact Bethe-Ansatz results (1D).0'0 
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{St = t) the exact W is recovered [Eq. (|TOD]. (ii) Three dimensional (3D) lattices are a 
further interesting direction for future developments. Indeed, encouraging results have been 
obtained for the simple cubic lattice at half-band filling. LDFT with Eq. ( |Tl| ) for W yields 
Egs/t = 1.21, 0.81, and 0.59 for U/t = A, 8 and 12, respectively, in good agreement with 
corresponding quantum Monte Carlo results,lii namely, E^^'^ /t = 1.27, 0.78, and 0.57. (iii) 
An accurate approximation to Wl'yu) has been also derived for the attractive (negative U) 
Hubbard model in an analogous way as for U > 0. For a ID ring with Na = Ng = 12 
we find \Egs - = 1.2 x lO^^, 7.5 x lO^^, and 1.5 x 10"^ for \U\/t = 1, 4, and 

64, respectively. These results show that LDFT describes electronic correlations correctly 
also when intra-atomic pairing is favored. Systematic investigations along these lines are 
currently in progress and will be published elsewhere.0 



V. CONCLUSION 



A new density-functional approach to lattice-fermion models has been developed that is 
by all means independent of the homogeneous electron gas. A simple approximation to the 
interaction-energy functional is derived for the Hubbard model, which provides with a unified 
description of correlations in all interaction regimes from weak to strong coupling. Results 
for the ground-state energy and charge-excitation gap of ID and 2D systems demonstrate the 
ability of lattice density functional theory to describe quantitatively the subtle competition 
between kinetic charge fluctuations and correlation-induced localization. The scope of DFT 
is thereby extended to the limit of strong electron correlations. 

Several interesting directions open up with potential implications in various related areas. 
For example, one may explore more general approximations to ^^[7] and one may apply 
the present approach to richer physical situations such as low-symmetry systems, disorder, 
magnetic impurities, or multiband Hamiltonians. These developments should be relevant 
to the study of lattice fermion models and also in view of a DFT description of strong 
correlations from first-principles. 
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